#pragma once

#include <iostream>
#include <fstream>
#include <iomanip>
#include "vect_3d.h"
#include "head.h"
// #include "class_ma.cpp"
//  #include "class_ma_Complex.cpp"
// #include "class_mac16.cpp"
#include "svd_double.h"
// #include "svdlongdouble.cpp"
// #include "svdcmpld.cpp"
#include "gaulen.h"
// #include "cgmom.cpp"
#include <string>
#include <chrono>
#include <ctime>
#include <limits>
#include <vector>
#include <cmath>
#include <Eigen/Dense>
#include <boost/math/special_functions/bessel.hpp>
#include <mkl.h>
#include <omp.h>

enum Grid_type
{
    StaggeredTartanGrids = 1,
    TartanGrids = 2,
    ChebyshevGrids,
    NewChebyshevGrids
};

enum RBF_TYPE
{
    IMQ = 1,
    Bessel = 2,
    Cubic,
    Quintic,
    TPS
};

struct ReUesData {
    std::vector<std::complex<double>> GG;
    std::vector<std::vector<double>> phi_mat_Pinv;
    std::vector<double> spx, spy, spz;
};


void printCurrentTime(string printinfo);

void generate_Gaussian_nodes();

double Cal_Fit_c(int Nums_1D, double length_lambda);

double factorial(int n);

void generate_tartan_interpolation_zeros(dvector* r, int N, double length,
    std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ);

void generate_tartan_interpolation_zeros(dvector* r, int Nx, int Ny, int Nz, double length_x, double length_y, double length_z,
    std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ);

void generate_Chebyshev_interpolation_zeros(dvector* r, int Nx, int Ny, int Nz, double length_x, double length_y, double length_z,
    std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ);

void generate_NewChebyshev_interpolation_zeros(dvector* r, int Nx, int Ny, int Nz, double length_x, double length_y, double length_z,
    std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ);

void generate_interpolation_zeros(dvector* r, int Nx, int Ny, int Nz, double length_x, double length_y, double length_z,
    std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ);

void generate_new_interpolation_zeros(dvector* r, int N, double length);

void generate_GG(int num_sp, double _k, Complex* GG, dvector* r, dvector rs);

void generate_the_Monomials_matrix_fortran_format(double* A, int N, double* spx, double* spy, double* spz,
    double c, double _k, double lambda, int RBF_type = IMQ);

void generate_the_Monomials_matrix_fortran_format(vector<double>& phi_mat, int N, double* spx, double* spy, double* spz,
    double c, double _k, double lambda, int RBF_type = IMQ);

void scattered_point_RBF_interpolation(Complex& Ga, double xx, double yy, double zz, Complex* GG, double* P[],
    double* a, int N, double c,
    double* spx, double* spy, double* spz, double _k, double lambda, int RBF_type = IMQ);

void scattered_point_RBF_interpolation(Complex& Ga, double xx, double yy, double zz, Complex* GG, vector<vector<double>> Phi_mat_Pinv,
    int N, double c,
    double* spx, double* spy, double* spz, double _k, double lambda, int RBF_type = IMQ);

double LagrangeBaseFunction(int idx, int k, double* tk, double t);

void scattered_point_lagrange_interpolation(Complex& Ga, double xx, double yy, double zz, Complex* GG, int numsp_1D,
    vector<double> zeroX, vector<double> zeroY, vector<double> zeroZ);

void scattered_point_lagrange_interpolation(Complex& Ga, double xx, double yy, double zz, Complex* GG,
    int Nums_1D_x, int Nums_1D_y, int Nums_1D_z,
    vector<double> zeroX, vector<double> zeroY, vector<double> zeroZ);

double new_test_scattered_lagrange(const double length_new, const double lambda, int grid_type, int Nums_1D, const double c_new,
    std::vector<std::pair<double, double>>& Length_OptimalC_list);

pair<double, double> New_test_scattered_RBF(const double length_x, const double length_y, const double length_z, const double lambda, int grid_type,
    int Nums_1D_x, int Nums_1D_y, int Nums_1D_z, const double c_new,
    vector<std::pair<vector<int>, vector<double>>>& Nums1Dxyz_Lengthxyz_OptimalC_list, int RBF_type = IMQ);

vector<vector<double>> Eigen_Cal_PseudoInverse(vector<double> phi_mat, int K);